Figure 5
4/12/24
Load data and libraries
##################
# LOAD LIBRARIES #
##################
library(tidyverse)
library(Seurat)
library(SeuratObject)
library(tidyseurat)
library(cowplot)
library(ggrepel)
#library(ggpubr)
#library(ggalluvial)
library(png)
library(grid)
library(patchwork)
#library(ggpuber)
library(openxlsx)
source("../../bin/spatial_visualization.R")
source("../../bin/plotting_functions.R")
#########
# PATHS #
#########
input_dir <- "../../results/08_spatial_dist/"
result_dir <- "../../results/09_figures/"
if( isFALSE(dir.exists(result_dir)) ) { dir.create(result_dir,recursive = TRUE) }
epi_clus <- "^5$|^6$|^7|^9" # non-filt
#############
# LOAD DATA #
#############
# DEGs_table <- read_csv(paste0(input_dir,"subset_100/DGEs_condition_wilcox_epi_100.csv")) %>% filter(groups == "DMPA")
DEGs_table <- read_csv(paste0("../../results/06_DGE_condition_st_data/","DGEs_condition_wilcox.csv"))
DATA <- readRDS(paste0(input_dir,"seuratObj_spatial_dist.RDS"))
img <- readPNG("../../resources/Schematic figure ST.png")Volcanoplot and IHC image (Figure 6A-C)
# width 170 mm, 6.7 inches
ord <- c("Superficial","Upper IM","Lower IM","Basal", "8","3","4","0","2","1")
col <- set_names(c("#F40014","#FF8A1E","#C77CFF","#A750B2"), ord[1:4])
###########################
# VOLCANO PLOT SUBMUCOSA #
###########################
ord <- c("Superficial","Upper IM","Lower IM","Basal")
DEGs_filt <- DEGs_table %>%
filter((grepl(paste0(ord, collapse = "|^"), .$layers))) %>%
mutate(layers = factor(.$layers, levels = ord)) %>%
filter(p_val < 0.099)
a <- col[DEGs_filt$layers]
set.seed(1);A <- DEGs_filt %>%
Volcano.fun_logFC(., "layers",
y.axis="p-value",
lab_size = 2.1, dot_size = .3,
up=c(.2, 0.05), down = c(-.2, 0.05)) + # labeling: (logFC, p-value)
ylab("avg. Log2 fold change") +
annotate("text", x = 3.9, y = .99, size = 4, col = "#CB251C", label = "Up in DMPA",fontface = "bold") +
annotate("text", x = 3.7, y = -.85, size = 4, col = "#3C71B5", label = "Down in DMPA", fontface = "bold") +
theme(legend.position = "none", #y,x
axis.text.x = element_text(colour = col, margin = margin(t = 10)),
#axis.ticks.x.bottom = element_line(colour = col),
plot.margin = unit(c(1, .5, .4, 1), "lines")) #t,r,b,l
A
#######################
# PLOTS VENN DIAGRAM #
######################
ord <- c("Superficial","Upper IM","Lower IM","Basal", "8","3","4","0","2","1")
sig_table <- DEGs_table %>%
filter(p_val_adj < 0.05) %>%
#filter((grepl(paste0(ord, collapse = "|^"), .$layers))) %>%
mutate(layers = factor(.$layers, levels = ord)) %>%
mutate(Direction = ifelse(avg_log2FC > 0, "UP", "DOWN"))
venn_input.fun <- function(sig_table, gr_col){
t <- sig_table %>%
# filter(grepl("_", .$layers)) %>%
group_by({{gr_col}}) %>%
nest() %>%
ungroup() %>%
arrange({{gr_col}}) %>%
mutate(layers = ifelse(grepl("\\d", .$layers), paste0("Clus_",.$layers), as.character(.$layers) )) %>%
mutate(genes = map( data, ~pluck(.x, "gene"))) %>%
mutate(genes = setNames(.[["genes"]], {{gr_col}})) %>%
mutate(direction = map( data, ~pluck(.x, "Direction"))) %>%
mutate(direction = setNames(.[["direction"]], {{gr_col}}))
return(t)}
t <- venn_input.fun(sig_table, layers) # cluster overlaps
# "./Figures/06/"
library(nVennR)
d <- plotVenn(t$genes[c("Basal","Lower IM","Superficial","Upper IM")],
sNames = names(t$genes[c("Basal","Lower IM","Superficial","Upper IM")]),
systemShow=F,
labelRegions=F, showNumbers=F,
showPlot =T,
opacity=0.1,
setColors = col,
outFile=paste0("./", "Venn_epi_sig.svg")
)
# d <- plotVenn(t$genes[-c(5:8)],
# sNames = names(t$genes[-c(5:8)]),
# outFile=paste0(params$fig.path, "Venn_epi_sig.svg")
# )
# showSVG(d)
knitr::include_graphics(paste0("./", "Venn_epi_sig.svg"))
#########################
# COMBINE PANEL A AND C #
#########################
# dev.new(width=6.7, height=3.2, noRStudioGD = TRUE)
A_B <- plot_grid(A, NULL, rel_widths = c(1,1), labels = c('A', 'B'), hjust = 0) # failed attempt at reading the svg file created by nVennR
a <- xml2::as_list( read_xml( paste0( "./", "Venn_epi_sig.svg")) )
l <- map_depth(a$svg$defs, 2,~pluck(.x, "path") ) %>% flatten()
l <- map(2:5,~a$svg$defs[[.x]]) %>% flatten() %>%
s <- get_shape(l[1] )
map_chr(l, ~attr(.x,"id"))
attr(a$svg$defs$symbol,"id")
id <- map_chr(map(2:5,~a$svg$defs[[.x]]), ~attr(.x,"id"))
id <- map_chr(l, ~attr(.x,"id"))
x <- attributes(l[[1]]$path)
id <- l %>%
map_chr(., ~attr(.x,"id")) %>%
set_names(seq_along(.), .)
id %>%
map(., ~get_shape(l[[.x]]$path) ) %>%
map(., ~as_tibble(.x, .name_repair="unique"))
annot_coord <- id %>%
map(., ~get_shape(a2[[1]]$svg[[.x]]) ) %>%
map(., ~as_tibble(.x, .name_repair="unique"))
myV3 <- createVennObj(nSets = 3, sSizes = c(rep(1, 8)))
myV3 <- plotVenn(nVennObj = myV3, systemShow=T, showPlot =F, labelRegions=F,
showNumbers=F, opacity=0.2, borderWidth=4,
outFile=paste0("./", "Venn_epi_sig.svg"))Plot dotplot and tissue
####################
# DOTPLOT FUNCTION #
####################
names <- c("Basal", "0", "1", "2", "3", "4", "8", "10")
clus_1 <- c("#F40014","#FF8A1E","#C77CFF","#A750B2") %>% set_names(., ord[1:4])
ID <- c("P118", "P105", "P080", "P031", "P097", "P108", "P114", "P107")
layer_dotplot.fun <- function(DATA, feat,
facet = TRUE,
txt_size = 5,
density = FALSE,
x_max=NULL,
col){
rects <- DATA %>%
group_by(layers) %>%
summarise(., ystart=min(spatial_dist, na.rm=T), yend=max(spatial_dist, na.rm=T),
Q1=quantile(spatial_dist, probs = 0.009, na.rm=T),
Q3=quantile(spatial_dist, probs = 0.90, na.rm=T)) %>%
filter(!(is.infinite(.$ystart))) %>%
mutate(Q1 = ifelse(.$Q1 == min(.$Q1), 0,.$Q1)) %>%
mutate(Q3 = ifelse(.$Q3 == max(.$Q3), max(.$yend),.$Q3)) %>%
mutate(Q3 = ifelse(.$layers == "Lower IM", .$Q3+.3,.$Q3)) %>%
#mutate(Q1 = ifelse(.$layers == "Basal", .$Q1+.3,.$Q1)) %>%
arrange(layers) %>% ungroup()
# rects <- DATA %>%
# group_by(layers) %>%
# summarise(., ystart=min(spatial_dist, na.rm=T), yend=max(spatial_dist, na.rm=T),
# Q1=quantile(spatial_dist, probs = 0.1, na.rm=T),
# Q3=quantile(spatial_dist, probs = .9, na.rm=T)) %>%
# filter(!(is.infinite(.$ystart))) %>%
# mutate(Q1 = ifelse(.$Q1 == min(.$Q1), 0,.$Q1)) %>%
# mutate(Q3 = ifelse(.$Q3 == max(.$Q3), max(.$yend),.$Q3)) %>%
# #mutate(Q3 = ifelse(.$layers == "Basal_2", .$Q3+.3,.$Q3)) %>%
# arrange(layers) %>% ungroup() %>%
# filter(., grepl(keep, .$layers))
DAT <- DATA %>%
#filter(., grepl("SubMuc", DATA$sp_annot) & !(grepl(epi_clus, DATA$Clusters))) %>%
#filter(., (grepl(keep, DATA$layers))) %>%
mutate(., FetchData(., vars = c(feat), slot = "data") ) %>%
select(groups, ID="orig.ident", layers, all_of(c(feat)), spatial_dist) #%>%
#mutate(ID = factor(.$ID, levels = ID))
dens <- ggplot() +
geom_density(data = DAT %>% filter(., .[[feat]] != 0 ), aes(y=spatial_dist, color=ID)) +
scale_y_continuous(expand = c(0, 0)) + ggpubr::clean_theme()+
scale_color_manual(values = col) +
theme(panel.border = element_blank(),
plot.margin = unit(c(0,0,0,0), "inches"))
if(facet == TRUE){facets <- facet_wrap(~groups, ncol = 2)
w <- 3}
else{facets <- NULL
w <- 1}
dot <- ggplot() +
#ggtitle(feature) +
geom_rect(data = rects, alpha = 0.1, show.legend=FALSE,
aes(xmin = -Inf, xmax = Inf, ymin = Q1, ymax = Q3, fill = layers)) +
geom_jitter(data = DAT, aes(x=.data[[feat]], y=spatial_dist, col=layers),
width = 0.1, alpha = 0.7, size=.3) + #
scale_fill_manual(values = col) +
scale_colour_manual(values = col) +
coord_cartesian(clip = "off") +
guides(col = guide_legend(override.aes = list(size=2), keyheight = .7, keywidth = .7, title = "Cluster")) +
scale_y_reverse(expand = c(0, 0)) +
#scale_y_continuous(expand = c(0, 0)) +
#scale_x_continuous(expand = c(0, 0)) +
{if(!(is.null(x_max))){xlim(-.5, x_max)}} +
facets +
my_theme + ylab("Similarity in gene expression") +
theme(legend.position="top", legend.justification = "left",
#legend.position = c(.1, 1.15), legend.direction = "horizontal",
plot.background = element_rect(fill = "transparent", colour = NA),
plot.margin = unit(c(2,.0,-2,.5), "lines"),
panel.border = element_blank(),
text = element_text(size = txt_size),
axis.line = element_line(linewidth = .3),
#axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.title.x.bottom = element_text(margin = margin(t = .5)),
axis.title.y.left = element_text(margin = margin(r = 5)),
panel.grid.major = element_line(linewidth = 0.2),
panel.grid.minor = element_line(linewidth = 0.1))
if(density == TRUE){dot <- dot + dens + plot_layout(ncol=2, nrow=1, widths=c(w, 1), guides = "collect")}
return(dot)
}###############
# FILTER DATA #
###############
DAT <- DATA %>%
filter(!(grepl("\\d", .$layers))) %>%
mutate(groups = if_else(groups =="ctrl", "Control", .$groups))
##############################
# CLUS DOTPLOT PER EPI LAYER #
##############################
# dev.new(height=1.9, width=6.7, noRStudioGD = TRUE)
# dev.new(height=2.4, width=6.7, noRStudioGD = TRUE) # with legend
feat <- c("IGHA1","SPRR2G","PKP1","MTM1" )
dot_fig <- map(feat,
~layer_dotplot.fun(DAT, .x,
facet = T,
txt_size = 9,
col = clus_1)) # , x_max = 6
# dot_fig[[1]]
C_1 <- wrap_plots(dot_fig, nrow = 1, guides = "collect" ) +
plot_layout(axis_titles = "collect") &
theme(#legend.position = c(.1, 1.15), legend.direction = "horizontal",
plot.margin = unit(c(1,.2,0,.3), "lines"),
legend.position = 'top',
legend.justification = "right",
legend.background = element_rect(colour ="gray", linewidth=.2),
legend.margin=margin(1,2,1,1),
legend.box.margin=margin(-25,0,-10,0),
#legend.box.margin=margin(-30,0,-15,0),
# axis.title.y.left = element_text(margin = margin(r = 5))
)
ggsave("./Figures/05/Fig_05C.png", C_1, width = 6.7, height = 2.4) #, dpi = 300
############################
# CLUS EXPRESION ON TISSUE #
############################
DAT <- DATA %>%
filter(grepl("P118|P097", .$orig.ident))
col_feat <- c("#EFEDF5", "#DADAEB", "#BCBDDC", "#9E9AC8", "#807DBA", "#6A51A3", "#54278F", "#3F007D") # Purples
p <- map(feat,
~plot_st_feat.fun( DAT,
geneid = .x,
txt = F,
save_space = F,
zoom = "zoom",
col = col_feat,
alpha = .9,
ncol = 2,
#annot_col = "#dbd9d9",
annot_line = .1,
img_alpha = 0,
point_size = .6))
# dev.new(height=1, width=6.7/4, noRStudioGD = TRUE)
# p_list[[1]]
# dev.new(width=6.7, height=1, noRStudioGD = TRUE)
p_list <- map(p, ~ .x + theme(plot.margin = unit(c(-25, -10, -10, -15), "pt"), #t,r,b,l
#panel.background = element_rect(fill = "white"),
panel.spacing.y = unit(-1.5, "lines"),
panel.spacing.x = unit(-2, "lines"),
legend.key.height = unit(0.2, "cm"),
legend.box.margin=margin(-20,0,0,0),
legend.title = element_blank(),
legend.position = "none"
) )
#########################
# COMBINE PANEL A AND C #
#########################
# patchwork does not respect margins when nesting plots, use cowplot instead!
# dev.new(width=6.7, height=3.4, noRStudioGD = TRUE)
c <- plot_grid( plotlist = p_list, ncol = 4, rel_heights = c(1), rel_widths = c(1))
C_2 <- plot_grid(NULL, c, rel_widths = c(.04,1))
(C <- plot_grid( C_1, C_2, ncol=1, rel_heights = c(3,1.1), labels = c('C'), hjust = 0))
# ggsave("./Figures/05/dotplot_gene_traject.png", C, width = 6.7, height = 3.2, bg = "white", dpi = 300)#############################
# COMBINE ALL FIGURE PANELS #
#############################
# dev.new(width=6.7, height=6.7, noRStudioGD = TRUE)
Figure5 <- plot_grid( A_B, C, ncol=1, rel_heights = c(1,1))
ggsave("./Figures/05/Figure05.pdf", Figure5, width = 6.7, height = 6.7, bg = "white", dpi = 1000)
Figure5